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DERIVATION OF NEWTONIAN-TYPE INTEGRATION COEFFICIENTS 
AND SOME APPLICATIONS TO ORBIT CALCULATIONS 


by 

C. E. Velez 
and 

J. L. Maury 

Goddard Space Flight Center 


CHAPTER I 
INTRODUCTION 


It is well known that finite difference operator techniques can be used to obtain many useful interpola- 
tion and integration formulas. In particular, recursive relations for the coefficients of various integration or 
quadrature methods can be obtained which are, in contrast to those obtainable by the method of undetermined 
coefficients, readily amenable to automatic computation. 

In this report, these difference operator techniques are used to construct generalized operators which 
define the coefficients of a large class of stable integration formulas of the Newton interpolatory type. The 
resulting recursive relations have been programmed and used to compute the coefficients associated with the 
various popular integration formulas such as those of Adams, Cowell, and Nystrom, as well as formulas of 
the Newton-Cotes type, which have applications in block, single, and multistep starting algorithms. A spe- 
cial application to the numerical integration of satellite orbits in multirevolution steps is also presented. 

Finally, a computer program which performs the calculations with rational arithmetic is described, and 
the coefficients associated with some of the well-known techniques are tabulated. 

It should be remarked that this report does not intend to present new formulas (although some of those 
derived are not easily found in the literature, in particular those pertaining to multirevolution starters) but 
to present a unified approach to many types of formulas which are currently being used to solve a variety of 
problems. 
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CHAPTER II 


FORMULATION 


A. Difference Operators and Identities 

Let n be a positive integer, s and h any real numbers, and f a real- valued function defined on an inter- 
val [a, b] that is n-times continuously differentiable; that is, f€C“[a, 5]. Consider the linear operators A^, 
E®, D", / defined by 


A„f(x) = fix + nh) - fix) 

(forward difference) , 

V„/(x) = /(X) - fix - nh) 

(backward difference) , 

Epix) ^ fix + snh) 

(shifting) , 

D"fix) = /(”>(x) 

(differentiation) , 

I fix) = fix) 

(identity) . 


and 


Interpreting equality between expressions containing these operators in the usual way,* some well 
known relations between these operators are the following: 

‘.-A- 


( 1 ) 


and 


From these identities, it follows that 


and 


= (/ + A) = £ . 

= [/ - = [/ - £-®] , 
hD = -log (/ - V), 

Using (4) and (5), we can immediately form the identity 






( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 

( 7 ) 


*Generally with respect to some class of polynomials. See Reference 1 concerning the calculus of finite differences. 


3 


where L is any positive integer and J and K are real. In this identity, powers of operators are defined in 
the usual way. 

Our goal is to find expressions for coefficients == y^(J, K, L) so that (7) can be expressed in powers 
of V as 




" 00-^ 

Yi(J. K, L)W 

. 1=0 




( 8 ) 


Once this is accomplished, we will show how (8) can be used to yield integration formulas for initial 
value problems of the form 

= /(x, y) , (9) 

with the initial values vix^), y'Cx^), . . . 

In similar fashion, using (3), (4), and (5), we have the identity 


L-l 





(Khpy 




CO 

r 

i~L 


(Khoy 




r 

L-iog(/- V)J 


E^hOy 


( 10 ) 


where L is any positive integer and J and K are real, and again we are to find expressions for coefficients 
= a^(J, K, L) so that (10) can be expressed as 


— - D^E~^ + aj(J, K, L)S/'D'^E^ 

i-l 1=0 


( 11 ) 


Equation (11) will be used to develop multistep starting formulas for systems of the form (9), especially for 
the case L - 2, which is of interest for orbit trajectory computations. 

Finally, for applications in the theory of multirevolution integration of satellite orbits, we require the 
operator identity 


V 


KN - 


(;-E/)g/ 
(' - '';") - ' 



( 12 ) 


which follows directly from relations (4) and (6). As before, we will seek an expansion of (12) of the form 




" °p ^ 

K, N)% 
i=0 




(13) 
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and will use it to develop multirevolution predictor-corrector and starting formulas. 

B. Series Expansions 

Before performing the expansions (8), (11), and (13), we will require the following series identities: 


where 


00 

[-log (1 - x)]^= , 

i^O 

1 i + 1 ’ 


for L > 1 , 






j=o 


i + L 


( 14 ) 


and 


00 

(l-x)”=y^(-l)^(”)x^ 

i=0 

l-(l-x)” = x^(-l)'(.yxV 


(15) 


(16) 


Identity (14) can be proved by induction as follows: For L = 1, 


[-log (1 


00 .• 

i=0 


is well known. Assuming the identity is true for L = n, we have 

[-log (1 - X)]”+1 = D-^D[-log (1 - x)]"+l , 

= 0-l)(n + l)[-log (1 - x)]"(l - x)-lj. 


By assumption, 


[-log (1 -x)]"+l = D“1 


(ji + 1)! x'+” y ^ //(”) 

i =0 ;=0 


( i + n + 1 ) 

i =0 ;=0 
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Hence, by definition of 


00 

[-log (1 - x)]"+^ =(n + l)!x”+^ y ^ 


and hence (14) is true for all L. Identity (15) is well known and (16) follows directly from (15). 

We now proceed to expand the operator expression (7) in powers of V. Using (16), we first have 


/ - E~^ = 7 - (7 - V)^ = 


^ = V ^ b .(7f)V' . 
i =0 


where 




SO that the factor jy - £ ^ can be expressed, by repeated series multiplication, as 

[/ - E~^Y = ’ 


i^O 


where, if L = 1 , 


and, if L > 1 , 


and 


e^(K, D-hy/Q, 


eo = r 


L k=0 L/pO 

where 7 ^ are dummy indices. Next, we have 


h 


72 


^h~ii 


^73-/2 ( ^^~ 7 l-i 


, 7 = 1.2, 3 




^7^ 


So by multiplying these last two factors, we get 


00 

(7 - E~^f E-J = ^ 7.(J, K. L)V' , 

i=0 


(17) 


(18) 


(19) 


( 20 ) 


( 21 ) 
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where 


(- iri(J_^ej{K. L ) . 

Finally, we have, from (14), 


[-log (/ - V)]^ = (22) 

2 = 0 


We can now form, 

by series division, the result 


where 

[- log (/ - V)] ^ ^ > 

2-0 

(23) 

and 

I'D = ) 

(24) 





and these are precisely the coefficients required in (8). In the following sections, this relation will be used 
to obtain various integration formulas. 


Next, we wish to expand identity (10) in powers of V. We first have that 

, IILJ ~ 

l-L 

and by (14) we have 


l-L l-L 


.(Khoy ^,,;[-iog(/-V)]* 

ir = ■ ” 7 ^ ’ 


[-log (I - V)]^ 

T\ 





SO that by collecting terms in powers of V, we have 


£ 

l=L 


(KhD)^ 


(^i(K, L)V^‘+^ , 
i=0 


(25) 
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where 


c^{K. L) = ^2 ■ 

;=o 

Next, we have, from (20), 

^-(J+K) _ ^ ^ (_ i)i('^+^)v' , 

i=0 

and hence the first factor on the right-hand side of (10) can be expressed as 


( 26 ) 


i-L 


(Khoy' 




_S£^ 

= 22 ’ 
i=0 


(27) 


where 


and 


(Jq - ^0’ 


i=o 


(28) 


Moreover, since 


V 


log (/ - V) 


/- E~^ T 

log (/ - V)J 


is just the expanded portion of identity (7) with /f = 1, J = 0, we have 


V 


log (/ - V) 


00 ^ 

y;yj(0,l.L)V^ 

i-0 


(29) 


where are given in (24). Hence, we can see that the Oj = a,(J, K, L), required in (11), are given by 


and 


“o “ ^oyo 

i 

<^i = y^/jYi-r 


(30) 


;=o 
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Finally, we wish to expand identity (12) in powers of V^. By (2), we can see that expanding 



E 


-J 

N 


in powers of is precisely the expansion of 


(/ - 


in powers of V so that, by (21), we have 


Finally, expanding the factor 


i=0 


i=l 


i=0 


where 




we obtain the required ° (5j{J , K, N) in (13), 


and 


/^o " ^0 ’ 

i 

)3,. = //J, K, l)-y^g/Af)|8j_j. 
/=1 


(31) 


Before proceeding to the applications of these operator expansions, we remark that the identities (7), (10), 
and (12), which we have expanded, were selected because of the specific applications we had in mind, 
otherwise their selection was arbitrary. Also, the methods used to obtain the expansions were essentially 
the same in all three formulas, and if the need arose for another formula derivable from Newtonian operator 
methods the same techniques would be applicable. Finally, these “techniques” involve little more than 
some elementary series algebra and result in formulas amenable to automatic computations. 
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C. Integration Formulas and Measures of Accuracy 

Applying the operators (8) or (11) to an arbitrary, sufficiently smooth function y(x), we see that the re- 
sulting relation expresses differences of values of the function in terms of differences of its Lth derivative. 
For example, by (8), we have 



K, L)Vy ^>(x + Jh ) . 


(32) 


This is precisely the type of relation one needs to numerically integrate initial value problems of the form 
(9). Note however that Equation (8) is valid only with respect to a class of polynomials, in which case the 
sums are finite, and that for an arbitrary function y(x), even when sufficiently smooth, the corresponding 
series (32) may fail to converge. These same remarks hold for Equation (11): hence we wish to find an 
expression for the error resulting from the truncation of (8) or (11) after n terms when applied to such a func- 
tion. To this end, we wish to estimate the difference operators (Reference 2) and G^, defined by 


and 


L^[y(x)] = V|y(x) - ^ y/J, K, L)V'y(^>(x + Jh) , 

i=0 

Gjy(x)] = V^y(x) “ ~ ^ ’ 

i-1 i=0 


(33) 


(34) 


where y(x) is assumed, for convenience, to be an infinitely differentiable function defined on some interval 
[a, b] with the property that for any x and x + rh, contained in [a, h], and n > 0, we have 

A"y(")(x +rh)=\2 • (35) 

m~n 


Next, by using the identity 


m 

V”y(x) = (- iy(^y(x - ink) , 


i=0 

where n, m are arbitrary positive integers, we see that (33) and (34) can be expressed, in ordinate form, as 

L 


^[y(x)l = (- D'(t)y(* - 

i=0 i=0 


(36) 
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and 


where 


y f^U\i — JL— 

G|,[y(x)] = y(x) - y(x - Kh) - ^ ^ - Kh) - ^ ^ ^/a)y^^Hx + (J - i)h] , 

i-i i-U 

^i(y) = (-l)'^ ('”)/„ 
ni-i 


Further » applying the Taylor expansion (35), we get 


and 


( 37 ) 




— (- iK)^ 


m\ 




X) 


n r no 

Xf<« E 

i=0 i/]=L 


(J - i)^~^ 
{m - L)T 


h™y(™)(x) 


m =0 




m! 


i-o 


00 


/n=L Li-0 


, ^^(y)(J - 
(m - L)! 


m~L 


/]^y(^)(x) . (38) 


Similarly, 


Gft[y(x)] r. y(x) - 2^ ^ _ .^, h’^y<-'”\x) 

m-Q 1 = 1 m=i 


^ r ~ 00 

i=0 m=L 


(j - 1 )'"“^ 
(m - L)! 


A'”y^'”)(x) , 


8^K’"h'”y^'”\x) 

m=l 



^/a)(J - i)”-^ 
(m - L)! 


!/]^y(^)(x) , 


(39) 


where 


d 


m 


m 

E 

i-n 


;!(m - ;)! ’ 


for 1 < ra < L - 1 
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I 



and 




T 77 for m > L - 1 . 

r!(m - ;)! 


So we see that the operators L^, and G|, can be expressed as a series in powers of og 


00 


00 

Gjy(x)] = ^ 


B h”'yOn\x ) , 


where 


(-!)'©(- 


for m < L - 1 


m , 


y!(m - j)\ 


^ > in\ / j (m " L)! 


6™ = - 


^ (- ^i(a)(J - 1)”-^ 

/ ^ ;Km - ;■)! ^ j (ra - L)! 


for m > L - 1 . 


Now the operators L^[y(x)] and G^[y(x)] are said to be of order p if 


Lft[y(x)] = 0(hP+^) 
G Jy(x)] = 0(hP+^) 


+ higher order terms , 


which will be the case if and only if 


= 0 • for 0 < m < p + L 


in the expansions (40). 





I 


The order of an integration formula of the form (32), with a finite number of terms, can be defined as 
the order of its associated difference operator. In the following, we will adopt this convention, and we will 
assume the truncation error associated with formulas derived from (8) or (11) in the form 

or 


[41(a)] 


where ^ is some point in the integration interval, so that we are essentially neglecting the higher order 
terms in the error expansion by assuming that some generalized mean-value theorem is applicable (see for 
example Reference 2, p, 247). 

In a similar fashion, we would like to measure the error resulting from the truncation of operator (13) 
after n terms when applied to a sufficiently smooth function. We remark that although this formula does not 
relate functional values to differences of a derivative of this function it could be considered an integration 
formula if one considers the first order approximation of hO as A. In actual practice, this formula will be 
used in conjunction with an independent integrator which computes the necessary functional values for the 
right-hand side of (13). The details of such an application are described in Chapter III, Section C of this 
report. At this point we simply define a measure of accuracy of such a formula in precisely the same fash- 
sions as was done for (8) or (11). Thus, we wish to estimate the difference operator 

n 

H^[y(x)] = V^^y(x) - N ^ K, A0Vj,Ay(x + JNh) , (42) 

i=0 

where again y(x) is assumed to be sufficiently smooth so that (35) holds. As before, we begin by expressing 
in ordinate form: 


n 

= y(x) - y(x - KNh) - N ^ ^,(|8)Ay[x + N(J - i)h ] . (43) 

1=0 


where 


^i(^)=(-i)'^(7K- 

m=i 


Expanding the terms as before, we get 
(- KN)^ 


Wi[y(x)l = 

jn=l 


m! 


A'”y('”)(x) - N ^ ^ ^ ^ 

i=0 m=l 


[N(J - i) + 1]"” - [N(J - i)]” 


m\ 


/j”^y(™)(x) . 
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So we see that can be expressed in powers of h^y^^\x) as 


H/,l>(x)] = ^ ^ , (44) 

ra=l 


where 




“ + AT [N(J 
i=0 


i) + 1]® - [N(J - i)Y 


The operator H^[y(x)] is defined to be order p if 

H^[y(x)] ~ + higher order terms . 

As for the integration formulas (8) or (11), the truncation error associated with formulas derived from (13) 
will be in the form 


where is the first nonzero coefficient in (44). 

We remark at this point that throughout this report, the error terms associated with various formulas for 
quadrature or integration derived from these operators will be omitted. For the methods given in the appen- 
dix, the error terms and orders presented were obtained directly from expansions (40) and (44), which are, in 
general, only estimates. No attempt was made to rigorously determine their sharpest form or estimate their 
magnitude. It was felt that such an analysis would, in general, be lengthy, difficult, and outside the main 
thoughts of this report which revolve about the idea of using computer-oriented arithmetic to derive useful 
numerical methods. Rigorous estimates can, of course, be found readily in the literature. 
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CHAPTER ill 


APPLICATIONS 


A. Multistep Quadrature and Integration Formulas 
for First- and Second-Order Systems 

We begin this section by indicating how the operator identity (8) can be used to define some well-known 
quadrature methods used to obtain approximations of integrals of the form 


Letting 



/(x) dx . 


(h-a) 


where n is some positive integer, and 


we seek the coefficients , so that 


X- a + i/i , 



/(x) dx - 


h 


i=0 




(45) 


To this end, let L = 1, J = 0 and K = n in (8) and, applying the operator to the function F(x), where 

F'(x) = f(x) . 

we have, retaining n terms and omitting the truncation error. 


V„F(x) = h ^ y^(0, n. l)V'F'(x) , 
i=0 


which can be rewritten for x = x^ as 


••u 


( 46 ) 
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We see that since (45) can be considered the ordinate form of (46), it is simply required that expressions for 
the )/| be obtained. From (24) we have 


and 


yo - « » 


;=1 


Now from (21) we have 

j=0 

= e^(n, 1) , 



and, by (14), 


Therefore 


//}!) 


1 

7+ 1‘ 


and 


ro = " 


7 i = ( 





and finally, converting (46) to ordinate form, we have 


H'i(n) = (-!)' ^(™)y;„(0,/J, 1) . 

m=i 


For example, for n = 2, we obtain 


)'0 = 2 , Yi 


1 

^2 = 3’ 


1 


W -±- \U - 3 . W -Jl 

0"a’ 


(47) 


(48) 
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and we have the well-known Simpson rule; 



f(x)dx =^^2 + + 


V- 


Similarly, for n = 3, we obtain 


yo-3. yi-~2’ 

w w 

”^ 0 - 8 ’ 18 ’ 


which yields Newton’s 3/8 rule: 


r^3 

Jxn 


^2=1- y3 = -|' 


^2 "8’ 3- 8’ 


f(x) dx ^ I/JU3 + 3/g + 3/^ + /q] , 


and in general, the coefficients defined by (47) and (48) yield the Newton-Cotes formulas of the closed type. 
The Newton-Cotes formulas of the open type, i.e., those which do not involve the ordinates at the ends of 
the integration interval, can be obtained from (8) by letting K = n, L - 1, J - - 1, and retaining n - 2 terms: 





f(x) dx - h 


y/- 1. n. 
i=0 


As before, we have from (24), 


and 


where 


/o - ^ 


y, = /,(-l.n. 1)-^ 
;=i 


All 

; + 1’ 


1 

” ^ V ’ since “ (“ 1)^”^ » L - I , 
;=o 
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and hence 




and 


yo 


Yi 


~ n 




J=0 7=1 

which, together with (48) with J = - 1, can be used to define formulas of the type 

n~2 


I f(x) dx = h ^2^ _ i_,- ) • 

•'xq i=0 


For example, for n = 3, we have 


Xo-3, yi - 


w w 

H'o-2’ ^1 2’ 


and 


r^3 

Jxo 


3/] 


/(x)dx - — + /j]. 


Likewise for n = 4, we get the formula 


f /(x)dx=^[2/3-/2 + 2/-J. 


(49) 


and so forth (see Reference 1, pp. 73-74). 

Next we consider some well-known integration methods derivable from (8) for initial value problems of 
form (9). For the case when L = 1, we wish to examine multistep methods of the Newton type in the form 


n 

Vj^y(x) = h ^ ^ y^(J, K, l)V-*y '(x f Jh). 

i=0 


First taking the values = we have 

Xo(-l,l,l) = l 


18 


and 


Now, 


but 


So we have 
and 


y,.(- 1,1. !) = /,•(- 1.1. D- 


Z-r 7 + 1 ■ 

7=1 


f,(- 1 . 1 . 1) - ^ (- !)'■“' (h)«, (I- 1) • 

J=0 

i 

= i>/l) . since ^ = (- 1)'“^' ; 

7=0 




ilifi = 0, 
/ 0 if 7 > 0 . 


yo(- 1,1.1) = 1 


y,(- 1.1.1) 


1 - 



lizL 

7 + 1 ’ 


which defines the well-known Adams- Bashforth integration predictor formula 


y(x) - y(x - h) = 7i ^ y^(- 1, 1. l)V'/(x - h) 
1=0 


( 50 ) 


( 51 ) 


for solving initial value problems 

y' = y) 

y(*o) = ^0 • 


( 52 ) 


In similar fashion, the often-used associated corrector formula, the Adams-Moulton, can be obtained with 
K" - 1, j = 0: we have 

yo(0, 1. 1) - 1 

and 

y,. (0.1.1)= f/O.l.l)-^^. 

}=l 


19 



where 




= e/l. 1) = b.(l) = 0. if j > 0. 


So we have 


yo(0, 1. 1) = 1 


YiiO. 1. 1) 


'^ZizL 
Au j + 1 


and we obtain the formula 


which is often used to solve (52). 


II 

y(x - h)~ h y/0, 1, l)V^/(x) , 


Another popular formula, known as Nystrom’s, can be obtained with K = 2, J 

yo(-l,2, 1) = 2 


Yii-hZ, 1) = fj(- 1,2,1) 


E liii 
j+1’ 


and, as before. 


but for any ; > 1, 


/,(- 1 , 2 > 1 ) = ^ (- iy->(^.}^ej{2, 1 ) , 
1=0 




We therefore have = 1 for all i > 1, and hence 


yo(-1.2, 1) = 2 


yj(-l,2,l) = l 


^ yj-j 

Z-t j + 1 
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which defines the formula 


n 

y(x) - y(x - 2h)~ h ^ ^ y^(- 1, 2, l)V'/(x - h) , 
i-0 


defining a predictor for (52). 

We now examine methods of the form (L = 2): 


n 

V|y(x) = ^ ^ YiiJ, K, 2)V'y"(x + Jti) , 


1=0 

for the solution of initial value problems of the form 

y"(x) = f(x. y) . 

y(Xo) = Yq . (55) 

y'(^o) = yo ’ 

which are of the type that frequently occurs in orbit trajectory computations. Analogous to the case when 
L = 1, we will obtain the predictor with J - - 1 and K - 1. We have, from (24), 


and 


Now 


yo(-1.1.2) = 1 

i 

y,.(- 1, 1, 2) = 1, 1.2)-^ 2!//f )y._. . 


7=1 


but 


7=0 

i 

7=0 


7 

e^.(l,2)=y^b^(l)bj_^(l) 

1=0 


2 ), 




/=o 
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jl if/ = 0, 
/0if;>0. 


Also 


so we see that 


Hf 


E 

l-n 


1 

I + 1 


i + 2 




'EM 


7i(-l,l,2) = l - 


,/=o 


7 + 2 


)^i- 


7=1 


are the coefficients defining the Stormer predictor formula 


for solving (55). 


n 

V^y(x) 

i-0 


y/-l,l,2)V^y"(x-/j) 


(56) 


Finally, in similar fashion, we can obtain the associated corrector formula, known as Cowell’s method, 
using J = 0 and K - 1: 


where 


and 


V^y(x) = yj.(0, 

i=0 


1. 2)V'y"(x) , 


yo(0.1.2) = l 


yi(0.1.2) = - 


I 



m=l 


(57) 


We remark here that recursions (50), (53), (54), (56), and (57) are well known (see for example Reference 
2)* and usually are obtained by other procedures. We note however, that the recursion (24) could be used to 
generate all these formulas in a unified fashion. 


*See also Maury, J. L., and Brodsky, G. P., '*Cowell Type Numerical Integration as Applied to Sotellite Orbit Computo- 
tion,” NASA X-553-69-46, April 1969. 
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We next proceed to other applications of the general operators (8), (11), and (13). The various integra- 
tion and quadrature formulas that can be generated by (8) through selection of values for J, K, and L should 
be clear at this point. It is noteworthy that by using nonintegral values of J and K, we can obtain a variety 
of interpolation methods which, together with the integration methods required here, could readily be used 
for stepsize modification or to compute output at nonstep points. 

It is also noteworthy that in addition to defining multistep methods, operator (8) can be used to define 
block-type Runge-Kutta methods (Reference 3); for example, the five-point formulas 

4 

Xi(- + Jh) , (58) 

i=0 

where K = 1, 2, 3, 4, and J = 4 - K, are precisely those required by such algorithms. 

Various combinations of the J , K, L parameters that were considered to be of general interest were 
used to compute the coefficients appearing in the appendix. 


B. Multistep Starting Formulas and Algorithms 

As is well known, a serious drawback of using multistep methods to integrate differential systems of 
the form (9) is the requirement that an independent method be employed to obtain the necessary starting val- 
ues. The problem is, in general, that the starting procedure used is not so efficient as the multistep inte- 
grator, requiring either more derivative evaluations per step or a smaller stepsize to maintain accuracy or 
both. In the present applications, in orbit and physical parameter estimation programs, this problem be- 
comes even more serious. This is so for the following reasons: 

(1) The force models employed in such programs are generally sophisticated, so that the running time 
of the program is directly proportional to the total number of derivative evaluations performed. 

(2) The number of equations to be integrated is frequently of the order of 100 or more. 

(3) It is often the case that the integrator must restart many times during the trajectory calculation be- 
cause of discontinuities introduced in the accelerations at discrete points in the orbit. These discontinui- 
ties may be due to thrusting, large solar radiation-shadow impulses, or a variety of other possibilities. 

The starting methods that are currently used in such applications are frequently based on either Runge- 
Kutta or power-series formulas. Another method that is known, but perhaps less frequently used, is based 
on multistep formulas. Algorithms based on these methods may often involve a more complicated computer 
program, but experimentation has indicated (Reference 4)* that such methods require fewer derivative evalu- 
ations than do the other methods. 

In this section, we present the formulas and describe the multistep starting algorithms for the case of 
second-order systems [L = 2 in (9)], since this case applies to orbit computations. The generalization of 
these ideas to more arbitrary differential systems should be clear from this example. Assume that we wish 
to apply the predictor -corrector formulas [see (56), (57)] 

*See also Peabody, P. R., “DODS Numerical Integration," Computer Sciences Corporation internal communication, 
September 1963. 
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(predictor) 


i=0 
k 

^^Va+I = y ] yi(0- 1’ (corrector) 

i=:0 


( 59 ) 


to integrate the differential system 


y" = f(x. y) 


with given initial values 


KXq) = 7o 


y'(Xo) = y o- 

Fixing k in (59), we can convert the formulas to ordinate form, 

i=0 

and 

k 

i=0 

We readily see that to start applying these formulas, we require the starting values 

|y,.,y"U = (0,k) (62) 

so that we could then use (60) to obtain y^^^ (n - k), followed by successive applications of (61) to obtain 
corrected values of y^^^ and thus completing the first step. We would then repeat this predictor-corrector 
cycle with N - k + 1, k + 2, and so on. Hence we are required to compute the values of (62) given Yq, y^, 
and yj = /(Xq, Yq). To this end we consider the operator (11), truncated after k terms, with L = 2, applied 
to y^; 

k 

i~0 

where the a^(J, K, 2) are given in (30). Letting K = 1, 2, 3, ... k md J - k - K in this expression, and con- 
verting each formula to ordinate form, we obtain the formulas 

k 

y^-yo = Oy()+ K=(l.k), (63) 

i=0 
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where, as before, 


^i(K) = (- 1)^' ^ (“)a^(fc - K, K, 2) . 
m=i 

The idea behind the method is to use (63) as one would a corrector formula, solving it by successive 
iterations. More precisely, let denote the mth approximation of the required values and assume 

these are known; then the (m + l)st approximation is given by 

k 

y^m+l) =y^ + Khy'Q + ^ ^ 

i=0 

and (64) 

, K - {l,k) , 

The convergence of such a scheme can be proved in the same manner as one does for the corrector formula 
(61) (see, for example. Reference 2, p. 216). In fact it can be shown that if C is the Lipschitz constant 
associated with /(x, y), or is bound on df/dy, then the method converges if h is such that 


max 

1< KSk 



< 1 . 


All that remains to be determined is a method to obtain the first approximation. In general, if h is suf- 
ficiently small, a sufficiently accurate first approximation is given by 


yP = yo+ Khy'Q+^^y'o) 

>/f = (l,fc). (65) 


In orbit trajectory calculations, a more accurate first approximation can be obtained by using an analytic 
two-body solution. These approximations are generally not expensive to compute and can considerably 
reduce the number of required successive iterations of (64). 

Another scheme which will generally reduce the number of required iterations of (64) is to place the 
given initial values in the center (assuming k is even) of the required starting values, so that they become 

( 66 ) 

The above algorithm can be used in this case with the changes K = (1, k) to K ~ (-k/2, k/2), K 0, and 
J=k~K to J = (k/2) - K, 


25 



We remark again that for higher values of L in (9), similar methods can be constructed from (11). Also, 
for L = 1, it is clear that (8) can be used to obtain starting formulas, in fact formulas (58) are precisely what 
we require for a five-point formula. Finally, we note that for equations containing derivatives, for example, 

= /(X, y, y', ... y(^~^)) , r > 1 , (67) 

the methods for different values of L can be used simultaneously to obtain the necessary starting values. 
This same remark applies, of course, to the basic integration of (67) using methods derivable from (8). 

C. Multirevolution Integration, Multirevolution 
Starting Formulas, and Algorithms 

To improve the efficiency of lifetime study and long-range prediction calculations, a method has been 
proposed that integrates orbits in multirevolution steps. This method is well known (for example, Reference 
5) and will be only outlined in this section. The recent interest in using orbit generation programs for life- 
time studies and planning interplanetary missions prompted the inclusion of these multirevolution integra- 
tion techniques in this report. Since starting procedures for such methods have not appeared in the popular 
literature, these are also included in the analysis. 

Essentially, the method of multirevolution integration involves combining a usual short-step numerical 
integrator with a procedure which steps the calculations ahead in multirevolution increments. This stepping 
procedure is similar to the usual predictor-corrector process in that it extrapolates the orbital elements N 
revolutions ahead and then, starting with these extrapolated values, computes successive corrections. 

We begin by outlining this basic algorithm (a more detailed description was given by Velez*). Let /y 
denote the value of an orbital element at the descending node of the yth revolution, and let N be the number 
of revolutions to be stepped. If k is the order of the highest difference to be retained in (13), we have for 
K = 1, J - - 1, the multirevolution predictor, applied to 

k 

0+JV = ^ ^ 1, 1. W)V^(A/p , (68) 

where we are using 
that is, h ~ 1 revolution, so that 


and where the are given by (31) as 

/So(-L-LN) = 


~ Vi 0’ 

£^1(A^,)=:A/,_^, 


*Velez, C. E., “Numerical Integration of Orbits in Multirevolution Steps," NASA X-S42-67-34 1, January 1967. 
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and 


1 

1, - 1. AO - /f(- 1. 1, 1) - ^ ■ 


i=i 

Now, as in (50), we have 

= alii, 

so we have, using the definition of gj(N), 


and 


1 . 1 . AO = 1 - ^ iy+1 (7+f 

7=1 


The associated multirevolution corrector can be obtained from (13) with K - 1, J = 0 applied to 




i=0 


(0,l,iV)V7^(Af.^^) 


( 69 ) 


where, since 


we have 


0 ( 0 . 1 . 1 ) 


1 if i = 0 
0 if 2 > 0 ’ 


and 


Bq(0,1,N) = 1 

p,(0, = (-!)<"■ 


y=i 


Further, we can express (68) and (69) in ordinate form: 


0+iV 


K 

2 = 0 


and 


K 

i=0 


(70) 


(71) 


27 



where and ^ are defined in the usual way in terms of the The basic algorithm can now be simply 
described as follows: 

(1) Compute the starting values 

A/., i = 0,iV,2iV kN, 

(2) Using the value j ~ kN, use (70) to predict the values of the orbital elements at the descending 
node of the [{k + l)iV]th revolution. 

(3) Using the short step integrator, and starting with these extrapolated values, integrate one revolu- 
tion to obtain their values at the [(ir + l)iV + l]th descending node and compute the difference 

(4) Using this difference and formula (71), correct the values of successively to convergence, 

repeating step (3) for each iteration. 

(5) Repeat steps (2) to (4) with } = (k + 1)N, (fe + 2)N and so on. 

It is easily seen that efficient starting procedures are essential to the overall efficiency of the above algo- 
rithm, both for the short step integrations in step 3 and for the starting values of step 1. These integrators 
require a knowledge of the orbital elements at the descending nodes of the first kN + 1 revolutions. In lieu 
of performing a short step integration over all these revolutions, the following method, analogous to those 
discussed in Sectiop B, could be used. 

Consider the operator (13), truncated after k terms applied to 

k 

^KN^KN - ^ ^ ] ^i(J, K, ’ 

i=0 

where we let = 1, 2, 3, . . . , k and J - k ~ K, Converting each such formula to the ordinate form, we have 

k 

- ^0 = ^ ^ K = (72) 

2 = 0 


where, of course, 

k 

m) = (- !)'■ y^(j);8^(fc -K,K,N). 

m~i 

And again, as in the case of using (63), the idea is to solve (72) by successive iterations. If we let 
and denote the mth approximation of these values and assume they are known, the (m + l)st approxi- 

mation is given by 

k 
i=0 

and 

K = (l,k). 
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where for each fC , is obtained by the short step integration of the elements to this node, using as 

initial values Note that the are just the values of the orbital elements at the node of the epoch 

revolution and are hence known. 

The convergence of this scheme can be proved in precisely the same manner as for (64). The first 
approximation of the starting values can readily be obtained from a two-body solution, and the idea of plac- 
ing the /q, A/q at the center of the required starting values, to improve convergence, can easily be formu- 
lated. 

We see that for each iteration, the method requires at most the short step integration of k revolutions, 
so that the overall efficiency of the method would generally be considerably improved over the starting pro- 
cedure based on a short step integration over kN 1 revolutions, especially for large N. 
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CHAPTER IV 


PROGRAM DESCRIPTION 

The series expansions of these generalized operators were programmed in Fortran IV on the Univac 
1108 using a rational arithmetic package to eliminate the deterioration that would have occurred using float- 
ing point. In the program, each rational number was represented by two contiguous, double-precision, 
floating-point words containing integral values. Input to the program consists of values for the variables J, 
K, and L. The output is tables of coefficients for the nonsummed* difference form, the nonsummed ordinate 
forms, and, where applicable, the summed ordinate forms. The nonsummed ordinate form is used by the pro- 
gram to determine the truncation error that defines the order of the method. The program consists of a main 
routine which determines the generalized operator to be used, a subroutine to create the ordinate forms, sev- 
eral subroutines to calculate the truncation error and order, an assembly language format routine, and the 
rational arithmetic package. 

The rational arithmetic package consists of— 

(1) GCD-a function using Euclid’s algorithm to compute the Greatest Common Divisor of two numbers 

[ap ag] = GCD >0, 

where GCD - 1 if a^ and a^ are 0 or if a ^ or a g is not integral. 

(2) ADD-a subroutine that performs rational addition defined by the following algorithm: 







[/3g, dg] 



[Ug, dg] 



(3) SUB -a subroutine that performs rational subtraction defined by 


,( fig) Ug 

dj dg d^ dg dg 


*For a discussion of the summed form of the integration formulas, see Reference 1, page 327, or Maury and Brodsky, 
mentioned earlier. 
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CHAPTER V 


SUMMARY AND CONCLUSIONS 


The recursive relations (24), (30), and (31) defining the coefficients of the operator identities (8), (11), 
and (13), respectively, were found to be an effective means of obtaining the coefficients of a broad spectrum 
of quadrature, interpolation, and integration formulas of the Newtonian type. The methods used to obtain 
these relations were elementary and allow considerable flexibility of the basic operator identities, easily 
yielding formulas of the Newtonian type not derivable directly from the operators considered in this report. 

Although not every possible application of these operators was discussed in detail (e.g. applications 
to interpolation) the applications presented should give the reader a sufficiently broad exposure to the capa- 
bilities, and also the limitations, of the difference operator technique. 

Finally, it is remarked that the coefficients defining the Adams-Cowell integration formulas (50), (53), 
(56), and (57) and the starters (63) are currently being used successfully in our orbit determination systems. 
Also, the multirevolution algorithm was tested and found to be an effective means of saving computer time 
for long arc calculations. The results of these tests were reported by Velez and mentioned earlier. It is 
expected that the multirevolution starters presented in this report will improve this process considerably. 


Goddard Space Flight Center 
National Aeronautics and Space Administration 
Greenbelt, Maryland. March 27, 1970 
311-07-21-01-51 
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APPENDIX 


NUMERICAL EXAMPLES 


Tables of Coefficients in Rational Form 

The following tables of coefficients are grouped according to usage. (For example, a table of predictor 
coefficients and a table of corrector coefficients form a group, the tables of coefficients for a multistep 
starter form another, and so on.) Preceding each group is a brief discussion and the pertinent operator equa- 
tion. Each table within a group is preceded by the value of J, the shift exponent; K, the block step varia- 
ble; and L, the derivative order. Each table is composed of two sets of coefficients: a set for use in the 
difference formulation and a set for use in the ordinate formulation. The local error constant (see Equations 
141(a)]) precedes each set of coefficients for the ordinate form formulas. 


Table Group I 

The two tables in this group are the Adams- Bashforth predictor coefficients and the Adams-Moulton cor- 
rector coefficients. These are used to solve the first-order differential equation 

y' = f(x, y ) . 


and, with Equation (8), can be expressed as 


y(x) - y(x - /j) = /J 


10 

^ ^ yj(J, 1, l)V'y'(x + Jh) , 
i=0 


where J = -- 1, 0, and by (24), 


and 


yo(J,l,D- 1, 


y,(J, 1. 1) = {,(j. 1. 1) - ^ ; + 


;=1 


f, -(-1,1,1) = !' 

// 0 , 1 , 1 ) = 0 ^ 


>i= 1,2,3 .... 
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The ordinate 11-point formulas are then given by 


II mil 


where 


10 

y(x) - y{x- h) = h ^/y)y'[x + (J - i)h ] , 
i=0 

10 

fi(y) = (-!)'■ 1.1). 

m=i 

J = -1 K = 1 L = 1 


Ordinate form, order =11, 
local error constant 4 111 223/17 418 240 

^0 = 2 132 509 567/479 001 600 
f , = - 2 067 948 781/1 19 750 400 
^2 = 1 572 737 587/31 933 440 
^3 = - 1 921 376 209/19 958 400 
= 3 539 798 831/26 611 200 
= -82 260 679/623 700 

= 2 492 064 913/26 611 200 
^7 = -186 080 291/3 991 680 
= 2 472 634 817/159 667 200 
= -52 841 941/17 107 200 

= 26 842 253/95 800 320 


J = 0 K = 1 L = 1 



Difference form 

local 

Ordinate form, order =11, 
error constant -4 671/788 480 

/o = 

1 

^0 

= 26 842 253/95 800 320 

/] = 

-1/2 


= 164 046 413/119 750 400 

/2 ^ 

-1/12 

^2 

= -296 725 183/159 667 200 

/3 = 

-1/24 

^3 

= 12 051 709/3 991 680 

/4 = 

- 19/720 

^4 

= -33 765 029/8 870 400 

ys - 

- 3/160 

^5 

2 227 571/623 700 

/6 ^ 

- 863/60 480 

^6 

= -21 677 723/8 870 400 

y? = 

-275/24 192 


= 23 643 791/19 958 400 

/8 = 

- 33 953/3 628 800 

^8 

= -12 318 413/31 933 440 

/9 = 

-8 183/1 036 800 

^9 

9 071 219/119 750 400 

/lo =-■ • 

- 3 250 433/479 001 600 

^10 

= -3 250 433/479 001 600 


Difference form 


/O 

/I 

/2 

/3 

/4 

/5 

/6 

/7 

/8 

/9 


1 

1/2 

5/12 

3/8 

251/720 

95/288 

19 087/60 480 
5 257/17 280 
1 070 017/3 628 800 
25 713/89 600 


y 10 = 26 842 253/95 800 320 
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Table Group II 


The two tables in this group are the Stormer predictor coefficients and the Cowell corrector coefficients. 
The formulas are used to solve the second-order differential equation 

and, with Equation (8), can be expressed as follows: 

10 

y(x) - 2y(x - b) ^ y(x - 2b) ~ b^ ^ ^ y/J, 1, 2)V^"(x + Jb) , 

where J = - 1, 0 and, by Equations (24), 

yo(J,l,2)= 1, 


1=0 


and 


1 

y/J, 1,2) = //J, 1. 2) - 2//}2)y._.(J, 1, 2) , 

7=1 

1) 

>2 = 1,2,3 


i) 


/,•(-!, 1,2) 

/jCO, 1,2) □ 0, 

The ordinate 11-point formulas are then given by 

y(x) - 2y(x - h) + y(x - 2h) = ^ ^ ^i(y)y"U + (J - i)h \ , 

with ^-{y) as previously defined in Table Group I. 


10 


1=0 



J =-l 

K 1 

L . 2 

Difference form 

Ordinate form, order - 11, 
local error constant 4 671/78 848 

ro - 

1 

^0 

- 263 465 639/159 667 200 

Xl = 

0 


= -296 725 183/79 833 600 

X2 - 

1/12 


= 1 742 930 263/159 667 200 

^3 -- 

1/12 

^3 

- 424 402 351/19 958 400 

Y4 = 

19/240 


- 2 337 301 223/79 833 600 

>'5 == 

3/40 


= -1 155 556 697/39 916 800 

ye = 

863/12 096 


- 1 637 523 683/79 833 600 

y? = 

275/4 032 

^7 

= -29 064 973/2 851 200 

yg = 33 

953/518 400 

^8 

539 999 083/159 667 200 

y, - 8 

183/129 600 

^9 

= -53 797 223/79 833 600 

y ,0 = 3 250 

433/53 222 400 

^10 

- 3 250 433/53 222 400 
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Ill 


J=0 K = 1 L = 2 


Difference form 


Ordinate form, order = 11, 
local error constant —24 377/13 305 600 


yo = 

1 

^0 

= 3 250 433/53 222 400 

y \ = 

-1 


= 3 1 24 027/3 193 344 

>'2 = 

1/12 

^2 

= -57 128 921/159 667 200 

^3 = 

0 


= 16 745 741/19 958 400 

Y 4 = 

- 1/240 

^4 

= -88 645 069/79 833 600 

>^5 = 

- 1/240 

^5 

= 42 375 577/39 916 800 

^6 = 

- 221/60 480 

^6 

= -2 342 533/3 193 344 

y ? = 

- 19/6 048 


= 7 139 837/19 958 400 

>^8 = 

-9 829/3 628 800 

^8 

= -18 674 153/159 667 200 

Y 9 = 

- 407/172 800 

^9 

= 1 838 819/79 833 600 

y\Q = - 

330 157/159 667 200 

^10 

= -330 157/159 667 200 


Table Group III 

The two tables in this group are the predictor formula and the corrector formula coefficients. The for- 
mulas may be used to solve the third-order differential equation 


With Equation (8), these can be expressed as follows: 


10 

V3y(x) = y.(J, 1 , 3)V'y'"(x + Jh ) , 

i=0 


where J = -1,0 and, by Equations (24). 

70 ^ 1 , 3 )= 1 . 

j 

Yi(J, 1, 3) = f/J. 1. 3) - ^ 1, 3) , 

and where the are given in Equation (21). The ordinate 11-point formulas are as defined in Table Groups 
I and II. 
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I 


Difference form 


yo 

= 

1 

/i 

= 

-1/2 

y2 

= 

0 

/3 

= 

0 

y4 

= 

1/240 

ys 

= 

1/160 

/6 

= 

221/30 240 

y? 

= 

95/12 096 

/8 

= 

9 829/1 209 600 

y9 

= 

2 849/345 600 

yio 

= 

330 157/39 916 800 


J - 0 K 1 

Difference form 

yo - 1 

y, = -3/2 
/2 = 1/2 

Xs = 0 

X4 = 1/240 

X5 = 1/480 

/6 = 1/945 

X7 = 11/20 160 

yg = 47/172 800 
Yg = 19/161 280 
= 430/15 966 720 


L = 3 

Ordinate form, order = 11, 
local error constant 24 277/2 956 800 

^0 = 6 275 141/11 404 800 

= 10 485 877/79 833 600 

^2 = 16 745 741/13 305 600 

^3 = -167 287/63 360 

^4 = 50 087 159/13 305 600 
^5 =-50 469 451/13 305 600 
= 4 522 591/1 663 200 

^7 = -3 020 741/2 217 600 
^8 = 2 419 061/5 322 240 

^9 = -7 261 259/79 833 600 
^10 = 330 157/39 916 800 

L ^ 3 

Ordinate form, order = 11, 
local error constant -61/2 280 960 

^0 = 330 157/39 916 800 

= 36 662 533/79 833 600 
^2 -- 15 601 049/26 611 200 

^3 = -100 921/950 400 

^4 = 5 935/66 528 

= -757 019/13 305 600 
= 124 909/4 435 200 

^7 = -3 251/316 800 

^8 = 34 189/13 305 600 

^9 = -6 271/15 966 720 

^,0 = 439/15 966 720 


Table Group IV 

The tables in this group are the coefficients of the starter formulas that can be used to form the starting 
values 


(yi.y,')> i = -5,-4, ... -1, 1,2 5 
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required by the formulas given in Table Group I to solve the first-order differential equation 


y =/(x,y) 


with the initial value 


KXq) = 7o • 


These formulas can be expressed in difference form as 

10 

y(Xg + Kh) - y(xQ) = h ^ y^(J, K, l)VV'(x^ + Jh) . 
i=0 

where J and K assume the values 


J 

10 9 8 7 

6 

4 

3 

2 

1 

0 


[-5 -4 -3 -2 

-1 

1 

2 

3 

4 

5 


The coefficients y^(J, K, 1) are given by Equation (24), and ordinate 11-point formulas are as defined in 
Table Groups I and II. 


J == 10 K =~5 L = 1 


Difference form 


Ordinate form, order = 11, 
local error constant -202 025/33 320 128 


yo = 

1 

fo 

114 985/19 160 064 

y, ■-= 

75/2 


= -320 875/4 790 016 

72 = 

-1 525/12 

^2 

311 375/912 384 

X3 = 

6 175/24 

^3 

= -838 375/798 336 

^4 = 

-49 775/144 

^4 

= 2 325 625/1 064 448 

ys = 

92 785/288 


= -89 035/24 948 

=- 

2 543 875/12 096 

^6 

= 2 306 375/1 064 448 

y? = 

84 375/896 

^7 

= - 2 793 625/798 336 


3 955 625/145 152 

^8 

= 2 996 375/6 386 688 

X9 = 

184 625/41 472 

^9 

= -8 183 125/4 790 016 

xio =-- 

5 256 4 25/19 160 064 

^10 

= -5 256 425/19 160 064 
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J =9 


K =-4 


L - 1 


Difference form 


/o = 

1 

y, = 

28/1 

yj = 

-260/3 

V3 - 

156/1 

V4 = 

-8 114/45 

/5 = 

1 250/9 


yj = -67 192/945 
= 21 808/945 
yg - -8 449/2 025 
/9 = 7/25 

/lO = 547/93 555 


J =8 K 

Difference form 


y, = 39/2 

r2 = -219/4 

yg = 693/8 

Y4 = -6 747/80 
yj = 8 253/160 

y^ = -43 021/2 240 
= 497/128 

yg =-12 881/44 800 
yj = -25/3 584 

y^o = -1 851/1 971 200 


Ordinate form, order = 11, 
local error constant 61/93 555 

^0 = -367/467 775 

= 4 099/467 775 

^2 = -1 387/31 185 

^3 = 2 996/22 275 

^4 = -13 462/51 975 
fs = 424/155 925 

^6 = -85 226/51 975 
^7 = - 14 972/31 185 
4 =-216 617/155 925 
=-158 327/467 775 
^10 = 547/93 555 


= -3 L = 1 

Ordinate form, order = 11, 
local error constant -827/3 942 400 

ifj = 127/394 240 

= - 1 857/492 800 

^2 = 40 433/1 971 200 

^3 = -17 247/246 400 
^4 = 34 983/197 120 

^5 = -5 149/7 700 

=-847 491/985 600 
^7 = -42 903/35 200 
^8 = -773 809/1 971 200 
^9 = 1 613/98 560 

^,0 = -1 851/1 971 200 
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J = 7 


K =~2 


L = 1 


Difference form 


Ordinate form, order = 11, 
local error constant 263/7 484 400 


ro = 

1 

^0 

= 

-263/7 484 400 

ri = 

12/1 


= 

263/748 440 

72 = 

-91/3 

^2 

= 

-3 439/2 494 800 

^3 = 

125/3 

^3 

= 

793/623 700 

74 = 

-2 999/90 

^4 

= 

7 213/415 800 

75 = 

688/45 

^5 


-62 389/155 925 

76 = 

-13 613/3 780 

^6 

= 

-102 569/83 160 

7l = 

41/140 

^7 

= 

-252 449/623 700 

75 = 

119/16 200 

^8 


8 249/356 400 

79 = 

13/14 175 

^9 


-9 707/3 742 200 

7\0 = 

251/1 496 880 

^10 

= 

251/1 496 880 


J =6 


K =-l 


L = 1 


Difference form 


Ordinate form, order =11, 
local error constant -14 797/191 600 6^ 


^0 = 

1 

^0 

= 

14 797/95 800 320 

7l = 

11/2 


= 

-32 309/17 107 200 

72 = 

-149/12 

^2 

= 

1 746 433/159 667 200 

^3 = 

117/8 

^3 

= 

-163 459/3 991 680 

74 = 

-6 731/720 

^4 


3 216 337/26 611 200 

/5 = 

4 277/1 440 

^5 


-379 571/623 700 

/6 = 

-19 087/60 480 


- - 

14 296 081/26 611 200 

7? = 

-275/24 192 

^7 


1 394 959/19 958 400 

78 = 

-7 297/3 628 800 


= 

-493 837/31 933 440 

79 = 

-7/12 800 

^9 


292 531/119 750 400 

/lO = 

-90 817/479 001 600 

^10 


-90 817/479 001 600 
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1 


J = 4 

K = 1 

L = 1 




Difference form 

Ordinate form, order = 11, 
local error constant -14 797/191 600 640 



ro 

1 

^0 

90 817/479 001 600 



y^ 

-9/2 


= -292 531/119 750 400 



y2 

95/12 

^2 

= 493 837/31 933 440 



^3 

= -161/24 

^3 

= -1 394 959/19 958 400 




= 1 901/720 

^4 

= 14 296 081/26 611 200 



Xs 

-95/288 


379 571/623 700 



>^6 

= -863/60 480 


= -3 216 337/26 611 200 



y? 

-13/4 480 

^7 

163 459/3 991 680 



ys 

= -3 233/3 628 800 

^8 

= -1 746 433/159 667 200 



y9 

= -2 497/7 257 600 

^9 

32 309/17 107 200 



y}o 

= -14 797/95 800 320 

^10 

= -14 797/95 800 320 




J = 3 

II 

to 

L = 1 



Difference form 

Ordinate form, order = 11, 
local error constant 263/7 484 400 



yo 

1 

^0 

= -251/1 496 880 



y\ 

= -8/1 


= 9 707/3 742 200 



y2 

= 37/3 

^2 

= -8 249/356 400 



y3 

= -9/1 

^3 

= 252 449/623 700 



y4 

= 269/90 

^4 

= 102 569/83 160 



ys 

= -14/45 

^5 

= 62 389/155 925 



y6 

= -37/3 780 

^6 

= -7 213/415 800 



yi 

= -1/756 

^7 

= -793/623 700 



yg 

= -23/113 400 

^8 

= 3 439/2 494 800 



X9 

0 

^9 

= -263/748 440 



XlO 

= 263/7 484 400 

^10 

263/7 484 400 
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L = 1 


J = 2 K = 3 


Difference form 


Ordinate form, order =11, 
local error constant —827/3 942 400 


/O 

= 

1 

^0 = 

1 851/1 971 200 

/l 

= 

-21/2 

= 

-1 613/98 560 

/2 


57/4 

^2 = 

773 809/1 971 200 

/3 


-75/8 

^3 = 

42 903/35 200 

/4 


237/80 

it 

847 491/985 600 

/5 


-51/160 

^5 = 

5 149/7 700 

/6 


-29/2 240 

It 

-34 983/197 120 

/7 

- 

-13/4 480 

^7 = 

17 247/246 400 

/8 


-7/6 400 


-40 433/1 971 200 

/9 

= 

-7/12 800 

^9 " 

1 857/492 800 

/]0 

= ■ 

-127/394 240 

o 

11 

-127/394 240 


J = 1 K =4 L = 1 


Difference form 


Ordinate form, order = 11, 
local error constant 61/93 555 


/o 

1 

/I 

= -12/1 

/2 

= 44/3 

/3 

= - 28/3 

>'4 

= 134/45 

/5 

= - 14/45 

/6 

= -8/945 

/7 

0 

/8 

= 13/14 175 

/9 

= 13/14 175 

/lO 

= 367/467 775 


= -547/93 555 

= 158 327/467 775 
^2 = 216 617/155 925 
^3 = 14 972/31 185 
^4 = 85 226/51 975 
= -424/155 925 

= 13 462/51 975 
^7 = -2 996/22 275 
fg = 1 387/31 185 

= -4 099/467 775 
^10 = 367/467 775 
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J =0 


K =5 


L = 1 


Difference form 


Ordinate form, order = 11, 
local error constant -202025/38 320 128 


yo = 

1 

= 

-25/2 

/2 = 

175/12 

/3 = 

-75/8 

/4 = 

425/144 

/5 = 

-95/288 

/6 = 

-275/12 096 

y? = 

-275/24 192 

^8 = 

- 175/20 736 

yy = 

-25/3 584 

/lo 

1 14 985/19 160 064 


^0 = 5 256 425/19 160 064 
^, = 8 183 125/4 790 016 
^2 =-2 996 375/6 386 688 
^3 = 2 793 625/798 336 
=-2 306 375/1 064 448 
^5 = 89 035/24 948 

= - 2 325 625/1 064 448 
^7 = 838 375/798 336 

^8 = -311 375/912 384 

^9 - 320 875/4 790 016 

^10 = - 114 985/19 160 064 


Table Group V 

The tables in this group are the coefficients of the starter formulas that can be used to form the starting 
values 

(y,. yp . i = -5. -4. ... -1,1.2. .... 5 

required by the formulas given in Table Group II to solve the second-order differential equation 

y" = f(x. y) 

with initial values 

7(Xo) = yo . 

and 

y'(xo) = yo • 

These formulas can be expressed in difference form as 

10 

y(Xo + Kh) - yfxg) = Khy' + ^ y/J. K, 2)VV"(Xk- + Jti ) . 

1 = 0 

where J and K assume the values shown in Table Group IV, the y^(J, K, 2) are given by Equation (24), and 
the ordinate 11-point formulas are as defined in Table Groups I and II . 
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Illlllllllllll■■■llllllllllllllllllllllllllll^ 


J = 10 


Difference form 


/O = 

25/2 

ri = 

-250/3 

72 = 

5 875/24 

73 = 

-30 125/72 

/4 ^ 

133 375/288 

/5 = 

-43 975/126 

/6 = 

4 404 125/24 192 

77 = 

-4 680 6 25/72 576 

78 = 

4 191 125/290 304 

79 = 

-8 183 125/4 790 016 

7i0 = 

202 025/3 483 648 


J = 9 

Difference form 


7o = 

8/1 

7i = 

-152/3 

72 = 

416/3 

73 = 

-9 664/45 

74 = 

1 864/9 

75 = 

-40 616/315 

76 = 

9 784/189 

7/ = 

-181 096/14 175 

78 = 

3 346/2 025 

79 = 

-1 094/18 711 

7l0 = 

-124/93 555 


K =-5 L = 2 

Ordinate form, order = 11, 
local error constant 1 918 325/1 162 377 216 

fg = -77 425/38 320 128 

.f, = 62 875/2 737 152 

^2 = -1 539 875/12 773 376 
^3 = 208 625/532 224 

^4 = -5 942 875/6 386 688 
^5 = 10 314 625/3 193 344 
= 22 426 625/6 386 688 
^7 = 5 650 375/1 596 672 
fg = 21 348 625/12 773 376 
^9 = 21 621 125/19 160 064 
^10 = 202 025/3 483 648 


K = -4 L = 2 

Ordinate form, order = 11, 
local error constant -1 556/70 945 875 

= -52/467 775 

= 758/467 775 

^2 = -356/31 185 

,^3 = 8 368/155 925 

^4 = -6 584/31 185 
^5 = 280 124/155 925 
= 532 184/155 925 
^7 = 2 704/1 485 

fg = 23 756/22 275 
^9 = 122/1 701 

^]0 = -124/93 555 
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J =8 


K =-3 


L =2 


Difference form 


Ordinate form, order = 11, 
local error constant 247 319/1 793 792 000 


yo = 

9/2 

^0 

= - 1 063/3 942 400 

y\ = 

-27/1 


6 511/1 971 200 

^2 = 

549/8 

^2 

= -10 833/563 200 

^3 = 

-3 831/40 

^3 

1 029/14 080 

74 = 

12 711/160 

^4 

= -88 827/394 240 

75 = 

-4 425/112 

^5 

= 280 821/197 120 

76 = 

50 319/4 480 

^6 

= 4 345 149/1 971 200 

77 =■ 

-35 451/22 400 

^7 

= 464 187/492 800 

Xs = 

225/3 584 

^8 

7 443/71 680 

79 = 

17/7 040 

^9 

-529/78 848 

7l0 = 

1 693/3 942 400 

fio 

1 693/3 942 400 


J =7 


K =-2 


L = 2 


Difference form 


Ordinate form, order =11, 
local error constant 1 303/21 021 000 


yo = 

2/1 

^0 

= 

-263/1 871 100 

>^1 = 

-34/3 


= 

263/149 688 

X2 = 

80/3 

^2 

= 

-131/12474 

73 = 

- 1 502/45 

^3 

= 

159/3 850 

74 = 

2 1 17/90 

^4 

= 

-41 543/311 850 

75 = 

-5 651/630 

^5 

= 

111 973/124 740 

y& = 

1 466/945 

^6 


35 932/31 185 

y? = 

-119/2 025 

^7 


263/5 670 

yo = 

-103/113 400 

^8 

= 

3 587/623 700 

79 = 

589/3 742 200 

^9 

= 

- 707/534 600 

y}o = 

109/935 550 

^10 

= 

109/935 550 
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J =6 


K =-1 


L =2 


Difference form 


Ordinate form, order =11, 
local error constant 5 512 813/145 297 152 000 


70 = 
y^ = 
Y2 = 

/3 = 
Yi = 

ys = 
re = 
y? = 
/8 = 
X9 = 
yio = 


1/2 
-8/3 
139/24 
-2 333/360 

5 539/1 440 
-2 713/2 520 

275/3 456 
8 563/1 814 400 

6 533/7 257 600 
30 577/119 750 400 
87 299/958 003 200 


= -14 797/191 6-00 640 

= 90 817/95 800 320 

^2 = -1 763 939/319 334 400 
^3 = 166 919/7 983 360 

=-10 111 819/159 667 200 
= 31 494 553/79 833 600 

= 14 797/82 944 

^7 = -60 917/1 900 800 

^8 = 466 157/63 866 880 

f 9 = -79 829/68 428 800 

^10 = 87 299/958 003 200 


J=4 K=1 L=2 

Ordinate form, order = 11, 

Difference form g^^^r constant -5 512 813/145 297 152 000 


YO = 

1/2 

^0 ■ 

87 299/958 003 200 

Y} = 

-7/3 

= 

-79 829/68 428 800 

Y2 = 

103/24 


466 157/63 866 880 

73 = 

- 1 387/360 

^3 = 

-60 917/1 900 800 

74 = 

475/288 

II 

14 797/82 944 

75 = 

- 1 231/5 040 

^5 = 

31 494 553/79 833 600 

/6 = 

- 199/24 192 

^6 = 

-10 111 819/159 667 200 

77 = 

- 409/259 200 

^7 = 

166 919/7 983 360 

78 = 

-3 391/7 257 600 

II 

-1 763 939/319 334 400 

79 = 

-263/1 496 880 

^9 - 

90 817/95 800 320 

7l0 = 

-14 797/191 600 640 

fio "" 

-14 797/191 600 640 



J = 3 


K =2 


L = 2 


Difference form 


Ordinate form, order = 11, 
local error constant -1 303/21 021 000 


yo = 

2/1 

^0 = 109/935 550 

7l = 

-26/3 

= -707/534 600 

X2 = 

44/3 

4^2 = 3 587/623 700 

ra = 

-538/45 

^3 = 263/5 670 

r4 = 

409/90 

^4 = 35 932/31 185 

/5 = 

-71/126 

^5 = 111 973/124 740 

^6 = 

- 19/945 

=-41 543/311 850 

77 = 

-52/14 175 

^7 = 159/3 850 

78 = 

-23/22 680 

^8 = -131/12 474 

79 = 

-263/748 440 

^9 = 263/149 688 

7l0 = 

-263/1 871 100 

^10 = -263/1 871 100 


J = 2 

K = 3 L = 2 

Difference form 

Ordinate form, order =11, 
local error constant -247 319/1 793 792 000 

7o - 

9/2 

^0 1 693/3 942 400 

7, 

-18/1 

= -529/78 848 

72 - 

225/8 

^2 = 7 443/71 680 

73 

- 849/40 

^3 = 464 187/492 800 

74 - 

1 203/160 

^4 = 4 345 149/1 971 200 

7s = 

- 123/140 

^5 - 280/821/197 120 

76 = 

-141/4 480 

= -88 827/394 240 

77 = 

-129/22 400 

^7 = 1 029/14 080 

78 = 

-21/12 800 

^8 = - 10 833/563 200 

79 = 

- 299/492 800 

^9 = 6 511/1 971 200 

7l0 = - 

1 063/3 942 400 

^,0 = -1 063/3 942 400 
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J=. 1 


K = 4 


L 2 


Difference form 


Ordinate form, order - 11, 
local error constant 1 556/70 945 875 


/o = 

8/1 

^0 

= -124/93 555 

/I = 

-88/3 


122/1 701 

X2 = 

128/3 

^2 

= 23 756/22 275 

/3 = 

- 1 376/45 

^3 

= 2 704/1 485 

74 = 

472/45 

^4 

= 532 184/155 925 

/5 = 

-376/315 

^5 

= 280 124/155 925 

/6 = 

-8/189 


= -6 584/31 185 

77 = 

-104/14 175 

^7 

= 8 368/155 925 

/8 = 

-26/14 175 


= -356/31 185 

/9 = 

-34/66 825 

^9 

758/467 775 

/lo- 

-52/467 775 

^10 

-52/467 775 


J = 0 K=5 L = 2 


Difference form 


Ordinate form, order = 11, 
local error constant -1 918 325/1 162 377 216 


Yo = 25/2 

= -125/3 

Y2 = 1 375/24 

yj = -2 875/72 

= 3 875/288 

y5 = - 1 525/1 008 
y^ = -1 375/24 192 

Yj = -125/10 368 

yg = - 1 375/290 304 
y, - -6 625/2 39 5 008 
yio --77 425/38 320 128 


^0 = 202 025/3 483 648 

= 21 621 125/19 160 064 
^2 - 21 348 625/12 773 376 
^3 = 5 650 375/1 596 672 
= 22 426 625/6 386 688 
^5 = 10 314 625/3 193 344 
= -5 942 875/6 386 688 
^7 = 208 625/532 224 

fg =-1 539875/12 773 376 
^9 = 62 875/2 737 152 

^,0 = -77 425/38 320 128 
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